# include "scene.h"
# include <stdio.h>
# include <math.h>
# include <stdlib.h>


static short r_conf_grid(struct grid_config *gc, struct grid *g);
static short th_conf_grid(struct grid_config *gc, struct grid *g);
static short z_conf_grid(struct grid_config *gc, struct grid *g);
static int nearest_int(double x);

extern short conf_grid(struct grid_config *gc, struct grid *g)
{
 short e;
#ifdef TRACE_FILE
 fprintf(g->trptr, "Entering into 'conf_grid()'...\n");
#endif

 e = r_conf_grid(gc, g);
 if(e >= 0) e = th_conf_grid(gc, g);
 if(e >= 0) e = z_conf_grid(gc, g);

#ifdef TRACE_FILE
 fprintf(g->trptr, "...'conf_grid()' ends.\n");
#endif

 return e;
}



static short r_conf_grid(struct grid_config *gc, struct grid *g)
{
 double temp;

#ifdef TRACE_FILE
 fprintf(g->trptr, "Entering into 'r_conf_grid()'...\n");
#endif

# ifdef r_UNIF_GRID
 fprintf(g->lptr, "Uniform grids are being generated in r-direction.\n");
 fprintf(g->lptr, "ERROR:: in 'r_conf_grid()': No definition for generation of uniform-grids in radial direction.\n");
 return -1;
# else

 fprintf(g->lptr, "Non-uniform grids are being generated in r-direction.\n");

 gc->nr_center = nearest_int(2*r_MID1/(dr_CENTER+dr_MID1));
 gc->r_cd_center = (dr_MID1-dr_CENTER)/(gc->nr_center-1);
 temp = 2*dr_CENTER/gc->r_cd_center - 1;
 g->nr_nozzle = nearest_int(0.5*temp*(sqrt(1+8*r_NOZZLE/(gc->r_cd_center*temp*temp))-1));
 gc->r_cd_center = 2*(r_NOZZLE/g->nr_nozzle - dr_CENTER)/(g->nr_nozzle-1);
 gc->dr_mid1 = dr_CENTER + (gc->nr_center-1)*gc->r_cd_center;
 gc->r_mid1 = 0.5*gc->nr_center*(dr_CENTER+gc->dr_mid1);

 gc->nr_mid12 = nearest_int(2*(r_MID2-gc->r_mid1)/(gc->dr_mid1+dr_MID2));
 gc->r_mid2 = gc->r_mid1 + 0.5*gc->nr_mid12*(gc->dr_mid1+dr_MID2);
 gc->r_cd_mid12 = (dr_MID2-gc->dr_mid1)/(gc->nr_mid12-1);

 gc->nr_base = nearest_int(2*(r_BASE-gc->r_mid2)/(dr_MID2+dr_BASE));
 gc->dr_base = 2*(r_BASE-gc->r_mid2)/gc->nr_base - dr_MID2;
 gc->r_cd_base = (gc->dr_base-dr_MID2)/(gc->nr_base-1);

 g->nr_gap = nearest_int(2*(r_CYL-r_BASE)/(gc->dr_base+dr_WALL));
 gc->dr_wall = 2*(r_CYL-r_BASE)/g->nr_gap - gc->dr_base;
 gc->r_cd_gap = (gc->dr_wall-gc->dr_base)/(g->nr_gap-1);

 g->nr_solid = gc->nr_center - g->nr_nozzle + gc->nr_mid12 + gc->nr_base;
 g->nr = g->nr_nozzle + g->nr_solid + g->nr_gap;

# endif

#ifdef TRACE_FILE
 fprintf(g->trptr, "...'r_conf_grid()' ends.\n");
#endif

 return 1;
}



static short th_conf_grid(struct grid_config *gc, struct grid *g)
{
 double temp;

#ifdef TRACE_FILE
 fprintf(g->trptr, "Entering into 'th_conf_grid()'...\n");
#endif

# ifdef th_UNIF_GRID
 fprintf(g->lptr, "Uniform grids are being generated in theta-direction.\n");
 g->nth = (int)((double)360/dTHETA); //Total number of cells in theta direction.


 if(g->nth%2 != 0) fprintf(g->lptr, "No. of cells in theta-direction must be even!!\nIncreasing it by one to nth = %d\n", ++(g->nth));

 temp = (double)360/(double)g->nth; //Re-defining "d_theta": in degrees.
 g->dth = temp*acos(0)/90; //"d_theta" in radians.

 if(nTH_COMP >= 2*NthG+1) g->nth = nTH_COMP; // No. of CELLs chosen for computation keeping in mind symmetry of the problem.
 else{
      fprintf(g->lptr, "ERROR:: in 'th_conf_grid()': Defined number of computational cells = %d. However, it should be more than '%d'. Redefine the parameter 'nTH_COMP' to anything more than '%d', re-compile, and rum the program again. Exiting.\n", nTH_COMP, 2*NthG, 2*NthG);
      exit(1);
     }

# else

 fprintf(g->lptr, "Non-uniform grids are being generated in theta-direction.\n");
 fprintf(g->lptr, "ERROR:: in 'th_conf_grid()': incomplete definition...\n");
 exit(1);
# endif

#ifdef TRACE_FILE
 fprintf(g->trptr, "...'th_conf_grid()' ends.\n");
#endif
 return 1;
}

 //[Modify both the conf_grid.c and grid_gen.c accordingly]

static short z_conf_grid(struct grid_config *gc, struct grid *g)
{
#ifdef TRACE_FILE
 fprintf(g->trptr, "Entering into 'z_conf_grid()'...\n");
#endif

# ifdef z_UNIF_GRID
 fprintf(g->lptr, "Uniform grids are being generated in z-direction.\n");

 g->nz = (int)(h_CYL/dZ);
 g->dz = h_CYL/g->nz;

# else

 fprintf(g->lptr, "Non-uniform grids are being generated in z-direction.\n");

 gc->nz_top = nearest_int(2*(z_TOP-z_MID2)/(dz_MID2+dz_TOP));
 gc->z_mid2 = z_TOP - 0.5*gc->nz_top*(dz_MID2+dz_TOP);
 gc->z_cd_top = (dz_TOP-dz_MID2)/(gc->nz_top - 1);

 gc->nz_mid12 = nearest_int(2*(gc->z_mid2-z_MID1)/(dz_MID1+dz_MID2));
 gc->z_mid1 = gc->z_mid2 - 0.5*gc->nz_mid12*(dz_MID1+dz_MID2);
 gc->z_cd_mid12 = (dz_MID2-dz_MID1)/(gc->nz_mid12-1);

 gc->nz_mid1 = nearest_int(2*(gc->z_mid1-nz_UNIF_BOT*dz_UNIF_BOT)/(dz_UNIF_BOT+dz_MID1));
 gc->dz_unif_bot = (gc->z_mid1-0.5*gc->nz_mid1*dz_MID1)/(nz_UNIF_BOT+0.5*gc->nz_mid1);
 gc->z_cd_bot = (dz_MID1-gc->dz_unif_bot)/(gc->nz_mid1-1);

 g->nz = nz_UNIF_BOT + gc->nz_mid1 + gc->nz_mid12 + gc->nz_top;
# endif

#ifdef TRACE_FILE
 fprintf(g->trptr, "...'z_conf_grid()' ends.\n");
#endif

 return 1;
}



static int nearest_int(double x)
   {
    short flag = 1;

    if(x < 0)
        {
         x = -x;
         flag = -1;
        }

    return flag*((x-(int)x < (int)x+1-x) ? (int)x : (int)x+1);
   }
